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Abstract. The random quantum Ashkin-Teller chain is studied numerically by means of time-dependent 
Density-Matrix Renormalization Group. The critical lines are estimated as the location of the peaks of the 
integrated autocorrelation times, computed from spin-spin and polarization-polarization autocorrelation 
functions. Disorder fluctuations of magnetization and polarization are observed to be maximum on these 
critical lines. Entanglement entropy leads to the same phase diagram, though with larger finite-size effects. 

The decay of spin-spin and polarization-polarization autocorrelation functions provides numerical evidence 
of the existence of a double Griffiths phase when taking into account finite-size effects. The two associated 
dynamical exponents z increase rapidly as the critical lines are approached, in agreement with the recent 
conjecture of a divergence at the two transitions in the thermodynamic limit. 

PACS. 05.30.Rt Quantum phase transitions - 05.70. Jk Gritical point phenomena - 05. lO.-a Computational 
methods in statistical physics and nonlinear dynamics 


1 Introduction 

Classical and quantum phase transitions are affected dif¬ 
ferently by the introduction of homogeneous disorder. In 
the former, it is well established that, when no frustration 
is induced, disorder is a relevant perturbation at a critical 
point when thermal fluctuations grow slower than disor¬ 
der ones inside the correlation volume. It follows that the 
critical behavior is unchanged when the specific heat ex¬ 
ponent a of the pure model is negative [T] . This criterion, 
due to Harris, has been extensively tested on classical toy 
models such as the 2D Ashkin-Teller model [5] or the 2D 
g-state Potts model EUD- In the latter, disorder is relevant 
for q > 2 and the new random fixed point depends on the 
number of states q. 

Quantum phase transitions, i.e. transitions driven by 
quantum fluctuations rather than thermal ones, involve 
new phenomena. First, randomness can never be consid¬ 
ered as homogeneous because time plays the role of an 
additional dimension. Therefore, in contrast to the classi¬ 
cal case, even when random couplings are homogeneously 
distributed on the lattice, they are always infinitely corre¬ 
lated in the time direction. Indeed, the random quantum 
Ising chain in a transverse field (RTFIM), for instance, 
is equivalent to the celebrated McCoy-Wu model, a clas¬ 
sical 2D Ising model with couplings that are randomly 
distributed in one direction but perfectly correlated in 
the second one [ 5115 ] . As a consequence, scale invariance 
is broken even after averaging over disorder. The random 


quantum fixed point is usually invariant under anisotropic 
scaling transformations. Correlation length ^ and auto¬ 
correlation time grow differently when approaching the 
random quantum critical point: 

6 - r, ( 1 ) 

where z is the dynamical exponent. In the RTFIM, or 
in any model whose critical behavior is described by the 
same fixed point, this dynamical exponent increases alge¬ 
braically when approaching the critical point and diverges 
at the critical point. 

Another feature of quantum phase transitions in pres¬ 
ence of disorder is the existence of Griffiths phases [7]. 
In the paramagnetic phase, there may exist large regions 
with a high concentration of strong couplings which can 
therefore order ferromagnetically earlier than the rest of 
the system. Even though the probability of such regions 
is exponentially small, they can cause a singular behavior 
of the free energy with respect to the magnetic field in a 
finite range of values of the quantum control parameter. 
The region of the paramagnetic phase where this phenom¬ 
ena occurs is called a disordered Griffiths phase. A sim¬ 
ilar phenomena takes place in the ferromagnetic phase. 
The singular behavior is due to regions of the system with 
a high concentration of weak bonds at their boundaries. 
They are therefore only weakly coupled to the rest of the 
system and can order independently [8] . Because the tun¬ 
neling time of these rare regions grows exponentially fast 
with their size, they have a drastic effect on the aver- 
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age autocorrelation functions of the system. Instead of the 
usual exponential decay, the latter displays an algebraic 
decay m 

( 2 ) 

involving the dynamical exponent z. In classical systems, 
Griffiths phases usually consist in essential singularities, 
too weak to be observed numerically, apart with some 
long-range correlated disorder [^[TO]. 

The quantum Ising chain in a transverse magnetic field 
has been, by far, the most studied system undergoing a 
quantum phase transition. The mapping of this model 
onto a lattice gas of free fermions allowed for exact cal¬ 
culations in the pure case m- In the presence of random 
couplings, exact results are sparse m but the mapping 
still allows for an efficient numerical estimate of static, 
as well as dynamic, quantum averages M- The critical 
behavior is governed by an unusual infinite-randomness 
fixed point (IRFP) which has been extensively studied 
using a real-space renormalization group approach, the 
Strong-Disorder Renormalization Group (SDRG), first in¬ 
troduced by Ma and Dasgupta and later extended to 
the RTFIM by Fisher [TBlII^fTB] . The strongest coupling, 
exchange interaction or transverse field, is decimated by 
projecting out the Hilbert space onto the ground state 
of this coupling. Other couplings are then treated using 
second-order perturbation theory. Nevertheless, the method 
is believed to become exact as the IRFP is approached 
because the probability distribution of random couplings 
becomes broader and broader and therefore, a strong cou¬ 
pling is always surrounded by weaker couplings that can 
be treated perturbatively. The dynamical exponent z was 
shown to diverge at the phase transition. The relation © 
is replaced by ^ ^ (In^t)^/’^ with ip = 1/2. Autocorrela¬ 
tion functions decay as m 


A{t) - (Intj-^"" (3) 

at the critical point, while correlation functions C{r) dis¬ 
play a more usual algebraic decay with the distance r. The 
Ma-Dasgupta renormalization group allows for the exact 
determination of the magnetization scaling dimension and 
the correlation length exponent naiiz]: 

2xa- = Wjv = 2 - ^ ^ v = llip = 2. (4) 

The approach has been applied numerically to higher di¬ 
mensions [2nl[2n[22] . The IRFP of the RTFIM is quite 
robust: in contrast to the classical case, the random quan¬ 
tum g-state Potts chain falls also into this universality 
class for any value of q [23l[24]- 

In this paper, a model with a richer phase diagram is 
considered. The quantum two-color Ashkin-Teller model 
can be seen as two coupled Ising chains in a transverse 


field. The Hamiltonian is [57] 

^ -E -E 

i i 

- ^ (5) 

i 

where and are two sets of Pauli matrices. The 
model possesses two Z 2 -symmetries, corresponding to the 
invariance of the Hamiltonian under the reversal of all 
spins Ui (or Ti) and of both Ui and r^. The breaking of 
these symmetries can be monitored using the two order 
parameters 


M = P = (6) 

i i 

referred to as magnetization and polarization. In the pure 
case, i.e. Ji = J, Ki = K, hi = h and gi = g, the phase 
diagram involves several critical lines, as the 2D classical 
Ashkin-Teller model [5Kl[5Hl[5Rll2Hj . When K < J, the two 
Z 2 symmetries are simultaneously broken and the Ashkin- 
Teller model undergoes a single second-order quantum 
phase transition with the control parameter 6 = J/h. The 
scaling dimensions of magnetization, polarization and en¬ 
ergy densities vary along the critical line m- 

1 TT 71- 

Xa = -z, X^r = -z--T, Xe = - --r 7) 

8 8arccos(—e) 2arccos(—e) 

for e = KfJ S [—l/\/2;I]. For K > J, i.e. e > I, the 
critical line splits into two lines, both belonging to the 
Ising universality class = 1/8, Xa-r = x„i =1/16 and 
Xe = 1). These lines separate the paramagnetic (M = 
P = 0) and Baxter {M, P ^ 0) phases from an intermedi¬ 
ate mixed phase {M = 0, P ^ 0). 


In the following, the random Ashkin-Teller chain is 
considered. The four couplings Ji, hi, Ki and gi are ran¬ 
dom variables, though not independent but constrained 
by the relation Q 



( 8 ) 


where e is a site-independent fixed parameter. This model 
was first studied numerically by means of Density-Matrix 
Renormalization Group (DMRG) in the weak-disorder regime 
e < 1 m- As in the pure model, the system undergoes a 
single quantum phase transition with the control param¬ 
eter _ _ 

^ = lnJ-lnh. (9) 


^ The case where Ki/Ji = ej and Qi/hi = th are different 
was considered in [31]. At the infinite-randomness fixed point, 
both quantities are renormalized to the same value, e* = 0 
in the weak-coupling regime {ej,eh < 1) and e —>■ -l-oo in the 
strong-coupling one. Without loss of generality, one can start 
with ej — eh- The more general case where ej and eh are 
random variables and are allowed to take values both above 
and below 1 was also considered in m and leads to a different 
critical behavior at the multicritical point. 
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SDRG shows that the inter-chain couplings Ki and gi 
are irrelevant on the critical line (5 = 0, i.e. the random 
Ashkin-Teller model behaves as two uncoupled random 
Ising chains. The critical behavior is therefore governed by 
the Fisher infinite-randomness fixed point with the criti¬ 
cal exponents 0. However, for finite disorder strength, a 
strong cross-over is observed numerically between the pure 
fixed point and this infinite-randomness fixed point. The 
regime e > 1 of the random Ashkin-Teller model was only 
studied more recently using SDRG [8T1I32] . The phase di¬ 
agram is qualitatively the same as the pure Ashkin-Teller 
model, in particular the two Ising lines still meets at a tri- 
critical point located at (5 = 0 and e = 1. When this point 
is approached by varying 6, the scaling dimensions © 
of the infinite-randomness Ising fixed point are recovered. 
However, when approaching this point along the half-line 
(5 = 0 and e > 1, the critical behavior is governed by 
different exponents: 


6- 2\/5 

TTTT’ 


i + Vf' 


( 10 ) 


Note that the ratio jS/v is unchanged, a property some¬ 
times referred to as weak universality. Between the two 
Ising lines in the regime e > 1, SDRG indicates the exis¬ 
tence of a double Griffiths phase: magnetization behaves 
as in the disordered Griffiths phase of the random Ising 
chain but polarization as in the ordered Griffiths phase. 


The critical behavior is expected to be unaffected by this 
choice. Indeed, the probability distributions of h and g, ini¬ 
tially delta peaks, will become broader and broader under 
renormalization so that the same IRFP will be eventually 
reached. This choice was made to minimize the number of 
disorder configurations. If L is the lattice size, the number 
of Ji couplings is L — 1 with open boundary conditions so 
the total number of disorder configurations is 2^“^. For 
small lattice sizes, up to L = 16, the average over disor¬ 
der can be performed exactly and the possibly disastrous 
consequences of an under-sampling of rare events can be 
avoided m- This strategy is motivated by the fact that 
we are mainly interested in Griffiths phases, where the 
dominant behavior is due to rare disorder configurations. 
The drawback is that a precise determination of critical 
exponents is more difficult, in contrast to |30| where the 
sampling was limited to 10,000 disorder configurations, al¬ 
lowing for larger lattice sizes up to L = 32. We also made 
additional calculations for a lattice size L = 20 but with 
an average over only 50,000 disorder configurations, ran¬ 
domly chosen among the 524,288 ones. As we will see, this 
under-sampling leads to observable deviations. 

For simplicity, we have moreover restricted ourselves 
to the case 

J 2 = 1/Ji In J, = 0. (13) 

and we have chosen a strong disorder by setting = 4 
and J 2 = 1/4. The quantum control parameter is now 


In the rest of the paper, new data of both regimes e < 1 
and e > 1 obtained by DMRG are presented and discussed. 
While only the critical point was considered in [30) . we 
are interested in the out-of-critical region of the phase di¬ 
agram and especially in the Griffiths phases when e > 1. 
In the first section, details about the implementation of 
the model and the parameters used for numerical com¬ 
putations are presented. In the second section, the phase 
boundaries are determined using integrated autocorrela¬ 
tion times, and the disorder fluctuations of magnetization 
and polarization. They are compared with the behavior of 
the entanglement entropy of one half of the lattice with the 
rest of the system. In the third section, the spin-spin and 
polarization-polarization autocorrelation functions are an¬ 
alyzed more carefully. In particular, we are interested in 
the algebraic decay © signaling the existence of a Grif¬ 
fiths phase. Finally, a conclusion follows. 


2 Numerical details 

We have considered a binary distribution of the intra¬ 
chain couplings Ji'. 

p{Ji) = ^ - Ji) + - J 2 )] ( 11 ) 

and homogeneous transverse fields h and g. Equation (© 
now reads 


5 = —Inh. (14) 


The model was studied using the time-dependent Density- 
Matrix Renormalization Group algorithm [Mll^l^ . A 
rough estimate of the ground state is first obtained with 
the so-called Infinite-Size DMRG algorithm. Because the 
couplings are inhomogeneous, the system was grown by 
adding single spins to one boundary rather than inserting 
them between the two blocks. After this initial Infinite- 
Size step, the accuracy of the ground state is improved 
by performing four sweeps of the Finite-Size algorithm. 
Since disorder fluctuations dominate at the IRFP, quan¬ 
tum fluctuations are expected to be much weaker than in 
the pure Ashkin-Teller model. For the latter, the expected 
critical exponents were recovered by keeping of the order 
of m = 192 states when truncating the Hilbert space of a 
left or right block in the DMRG algorithm. For the ran¬ 
dom Ashkin-Teller model, we fixed the upper limit of this 
parameter to m = 64. The actual number of states was 
determined dynamically by imposing a maximal trunca¬ 
tion error: 10“^ during the initial Infinite-Size step, 10“®, 
10“^, 10“®, and 10“® during the four Finite-Size sweeps. 
Using these parameters, we were able to make calculations 
for a large number of quantum control parameters 5 for 
lattice sizes up to L = 20 for e > 1. Unfortunately, the 
Arpack library, used to determine the ground-state in the 
truncated Hilbert space, sometimes failed for some partic¬ 
ular disorder configurations. In these cases, the point, and 
not simply this disorder configuration, is discarded. For 
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e < 1, many calculations failed for L = 16. Only 27 values 
of the control quantum parameter could be completed for 
e = 1, mostly far from the critical point. Moreover, when 
successful, the computation takes a time which increases 
very fast for e < 1. Since the two Ising chains are uncou¬ 
pled at the fixed point, the Hilbert space becomes closer to 
a tensor product of the spaces of two Ising chains. There¬ 
fore, the number of states to be kept during the truncation 
process of the DMRG algorithm should be of the order of 
the square of the number of states necessary for a single 
Ising chain. For this reason, the largest lattice size consid¬ 
ered in the regime e < 1 is only L = 12. 

Average magnetization and polarization densities 

(m) = (0|cr2/2l0)^ (t) = (0^2/2^1/210): ( 15 ) 


3.1 Integrated autocorrelation time 

One of the properties that define criticality is that any 
characteristic length or time disappears at a second-order 
phase transition. Out-of-criticality, the exponential decay 
of average spatial correlation functions C(r) and autocor¬ 
relation functions A{t) provides respectively a correlation 
length ^ and an autocorrelation time . In a pure system, 
both quantities are expected to diverge as a critical point 
is approached. In the random case, a divergence of ^ and 
is expected in the whole Griffiths phase. However, in a fi¬ 
nite system, these divergences are smoothed and replaced 
by a finite peak. At large time t, connected autocorrelation 
functions A{t) are dominated by an exponential decay of 
the variable t/^f Consequently, their integrals behave as 


were measured at the center of the chain. |0) denotes the 
ground state and the over-line bar stands for the average 
over disorder. In order to measure non-vanishing averages, 
longitudinal magnetic and electric fields were coupled to 
the two boundary spins of the chain with the Hamiltonian 

Hi = Bal + Ealrl + Bal + Ealrl (16) 

to break the two Z 2 symmetries. The convergence of the 
DMRG algorithm is also faster when such boundary fields 
are imposed. Spin-spin and polarization-polarization con¬ 
nected autocorrelation functions, defined as 


Aa{t) = (0|cr£/2(^)'^L/2(0)|0) - ("^)^ (17) 

Acrrit) = (0|cr2/2(l)^L/2(^)'^L/2(0)'^2/2(0)|0) “ (T’)^ 

were estimated using a discretized imaginary-time evolu¬ 
tion operator: 


Aa{nAt) 


(0|a£/2(l-gZ\t)"a£^2l0) 


{my. (18) 


We have used the value At = 10 ^ and computed auto¬ 
correlation functions up to t = 10. 


3 Phase boundaries 

As discussed in the introduction, the random quantum 
Ashkin-Teller model is expected to undergo a single tran¬ 
sition when e < 1 and two transitions when e > 1. This 
is easily observed on the behavior of magnetization and 
polarization, which are the two order parameters of these 
two transitions. As seen on figures [TJ magnetization and 
polarization display a fast variation but at different val¬ 
ues of the transverse field h, and therefore of the control 
parameter <5 = — In/i, when e > 1. 

However, because of the hnite-size of the system, mag¬ 
netization and polarization curves are too smooth to pro¬ 
vide accurate estimates of the location of the transitions. 
Diverging quantities are more convenient for that purpose 
and usually preferred in numerical studies. In this section, 
we discuss three quantities that diverge, or display a pro¬ 
nounced peak, at the transitions of the random Ashkin- 
Teller model. 


r+oc _ ^-1-00_ 

T = / A{t/^t)dt = / A{u)du (19) 

JQ Jo 

and, like should display a peak. We have computed 
the integrated autocorrelation time r for spin-spin and 
polarization-polarization autocorrelation functions. The up¬ 
per bound of the integral (fTOl) was replaced by the largest 
time t = 10 considered. This approximation has no effect 
on the estimate of the autocorrelation time r as long as 
is much smaller than 10. As will be seen below, this is 
the case for the lattice sizes that we considered. 

As can be seen on figures [5] and O the integrated au¬ 
tocorrelation times display two peaks. The first peak oc¬ 
curs at a value of the transverse field h which is of the 
same order of magnitude as J 2 . Therefore, this peak is 
probably associated to the ordering transition of the dis¬ 
order configurations with a majority of weak couplings 
J 2 . However, the height of this peak does not increase 
significantly with the lattice size so one can conjecture 
that this peak will remain finite in the thermodynamic 
limit and is not associated to any phase transition. The 
height of the second peak clearly increases with the lat¬ 
tice size. For e < 1, the location of the peak is roughly 
the same for spin-spin and polarization-polarization au¬ 
tocorrelation times. For e > 1, the data clearly shows 
that the peak occurs at a positive control parameter 5, 
i.e. a transverse field h < 1, for spin-spin autocorrela¬ 
tion functions and negative for polarization-polarization 
ones. This indicates that the system undergoes an electric 
phase transition followed, at larger control parameter, by 
a magnetic one. This is consistent with the picture given 
by magnetization and polarization curves. The location of 
the two transitions was predicted by Hrahsheh et al. to be 
(5c = ± In f for e » 1 m- For e = 4, we observe the two 
peaks at Sc = — In he — 0.54 and 6c — —0.99 for L = 16 
for instance, still far from ±ln| ~ ±0.69. Moreover, our 
transitions lines are not symmetric with respect to the axis 
(5 = 0, as required by self-duality. Since the data was pro¬ 
duced by DMRG with a relatively large number of states 
and since the averages were made over all disorder con¬ 
figurations, the deviation can only be the consequence of 
the relatively small lattice sizes that could be reached and 
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Fig. 1. Magnetization (ieft) and poiarization (right) of the random quantum Ashkin-Teller chain versus the transverse field h. 
The different curves correspond to different values ai t = KijJi. The lattice size is L = 12. 
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Fig. 2. Autocorrelation time estimated by integration of the average spin-spin autocorrelation function Acr{t). The different 
graphs correspond to different values of £ and the different curves to different lattice sizes L. 
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Fig. 3. Autocorrelation time estimated by integration of the average polarization-polarization autocorrelation function 
Ac,T{t). The different graphs correspond to different values of £ and the different curves to different lattice sizes L. 


of the boundary magnetic and electric fields which favor a 
Baxter phase and therefore shift the whole phase diagram. 

We also considered the first moment 

£- 1-00 _ £- 1-00 _ 

/ tA{t)dt / / A(t)dt (20) 

do do 

that should be equal to the autocorrelation time if the 
connected autocorrelation function A(t) displays a purely 
exponential decay A{t) ^ Like the autocorrela¬ 

tion time, the first moment was computed for both spin- 
spin and polarization-polarization autocorrelation func¬ 
tions. When plotted with respect to the transverse fields, 
two peaks are observed. Even though the shape of these 
peaks is not strictly identical to that of the autocorre¬ 
lation time (HU), in particular the second peak is higher 
and slightly broader, both quantities behave in the same 
way with the transverse field h. Therefore, the same con¬ 


clusions can be drawn. A reconstructed phase diagram is 
shown on figure|H It is qualitatively similar to the one pre¬ 
sented in Ref. |31) . However, it is not symmetric under the 
transformation 5 O —5. As discussed above, finite-size ef¬ 
fects are here strengthened by the boundary magnetic and 
electric fields that globally shift the phase diagram. 


3.2 Disorder fluctuations 


In a random system, any thermodynamic average (A) is 
the result of a quantum average 

(A) = {ilj^[Ji,K,]\X\^o[JuK,]) 


( 21 ) 
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Fig. 4. Phase diagram in the parameter space (e, h) obtained 
from spin-spin (continuous lines) and polarization-polarization 
(dashed lines) first moment (1201) . The different curves corre¬ 
spond to different lattice sizes. 

followed by an average over coupling configurations 
JX) = ^{1Po[J^,K,]\X\MJ^,K^])p{{J,,K,})Y[dJ,dK, 

i 

( 22 ) 

where Ki]) is the ground state of the system for 

a given coupling configuration {Ji,Ki} and p({Ji, Ki}) 
the probability of this configuration. At an IRFP, disor¬ 
der fluctuations dominate over quantum fluctuations. The 
strength of the former can be measured by the variance 

Vx=JW-V0"- (23) 

We computed this quantity for both magnetization (V)r) 
and polarization (V)rr)- As can be seen on figures [S] and 
[6l the variances 14 and Var are numerically very stable. 
They vanish at high and low transverse fields h and dis¬ 
play a well-defined single peak. In particular, there is no 
second peak at ^ J 2 . The locations of the maxima of 
the peaks are accurately determined and are in agreement 
with the ones estimated from autocorrelation times. The 
same conclusions can be drawn: the magnetic and electric 
transitions occur at very close control parameters 6, prob¬ 
ably the same, for e < 1, while a finite shift is observed 
for e > 1. Even though only a weak dependence on the 
lattice size L of 14 and t4r is observed on figures [5] and 
1 a systematic finite-size shift is present. For e < 1, the 
distance between the two critical lines decreases when the 
lattice size L increases, in agreement with the prediction 
of a unique transition. The coincidence of the maxima of 
the autocorrelation times with those of the disorder fluc¬ 
tuations shows that the phase transition is induced by 
disorder fluctuations, rather than quantum fluctuations, 
as expected at an IRFP. 

As can be noticed on figures [5] and [6l the height of the 
peaks of the variance of disorder fluctuations increases 
slightly with the lattice size, at least for L < 16. The 


data for the lattice size L = 20 display indeed a smaller 
peak. This lattice size is the only one for which the av¬ 
erage has not been computed over all possible disorder 
configurations but only over a subset (^ 10%) of them. 
The smaller peak for L = 20 is therefore probably due to 
an under-sampling of the dominant configurations at the 
critical point. 50, 000 is still, at least for certain quanti¬ 
ties, a too small number of disorder configurations. In the 
following, data for L = 20 should be taken with more care 
than smaller lattice sizes, for which an exact average over 
disorder was performed. 

3.3 Entanglement entropy 

When the degrees of freedom of the system can be divided 
into two subsets A and B, and therefore when the Hilbert 
space can be decomposed as a tensor product H = Ha <8 
Hb, the degree of entanglement of the two sub-blocks is 
conveniently measured by the von Neumann entanglement 
entropy of A with the rest of the system m- 

Sa =-Tr-HAPA'^ogpA (24) 

where pA is the reduced density matrix 

Pa = T:t:-HbP (25) 

and p the density matrix of the full system. In the case of 
a pure state IV’), the latter is the projector p= In 

the following, we will consider the subset A made of the i 
spins at the left of the chain. 

Entanglement entropy has recently attracted a lot of 
attention because of Conformal Field Theory (CFT) pre¬ 
dictions at pure critical points [55115^ . The predicted log¬ 
arithmic behavior with £ is also observed in RTFIM but 
with a prefactor that involves an effective central charge 
c = ^ In 2 [40] . The entanglement entropy is also com¬ 
monly used in the literature to determine phase bound¬ 
aries m- Indeed, it is expected to be larger when quan¬ 
tum correlation functions are long-range. At an IRFP, the 
entanglement entropy is related to the probability of a 
strongly correlated cluster across the boundary between 
the two blocks A and B. Numerically, the reduced den¬ 
sity matrix pA being computed and diagonalized at each 
step of the DMRG algorithm, the entanglement entropy 
is given without any additional computational effort. 


The average entanglement entropy S{£) of the random 
quantum Ashkin-Teller chain is plotted on figure |7| for 
I = L/2. For e = 4, two peaks can be observed and inter¬ 
preted as the signature of the two phase transitions. As 
expected, one single peak is present when e < I. However, 
only one peak can be distinguished in the case e = 2, while 
autocorrelation time indicates the existence of two transi¬ 
tions. Because of the finite-size of the system, the expected 
two peaks are probably merged into a single one. This sce¬ 
nario is compatible with what is observed for e = 4: what 
was only a shouldering at the left of the peak for L = 8 
becomes a second independent peak at L = 20. 
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Fig. 5. Variance of disorder fluctuations of magnetization. The different graphs correspond to different values of e and the 
different curves to different lattice sizes L. 


The phase diagram is qualitatively the same as previ¬ 
ously constructed. However, the two peaks are not located 
at the same position as those displayed by the integrated 
autocorrelation times, or the variance of disorder fluctua¬ 
tions. At e = 4, they are instead found at Sc ^ —0.10 and 
Sc ^ —1.33, far from the estimates Sc — 0.54 and —0.99. 
This large difference is probably due to Finite-Size effects. 
Indeed, magnetization, polarization and autocorrelation 
functions were computed at the center of the lattice, i.e. 
at the site L/2. In contrast, the entanglement entropy is 
a global quantity, therefore more sensitive to the presence 
of boundary fields. 

CFT predicts that the entanglement entropy of a block 
of size £ behaves as m 


S{£) = pc In 


■ L , Tr£' 
— sm — 
-7ra L. 


Cst., 


where c is the central charge and p is equal to 1/3 for pe¬ 
riodic boundary conditions and 1/6 for open boundaries. 
However, this relation was obtained on a finite but contin¬ 
uous manifold and not on a lattice. Therefore, it is only 
poorly verified by our numerical data, for which strong 
lattice effects are still present. Nevertheless, the predicted 
dependence on the lattice size L is well reproduced by the 
numerical data. For an equal partition of the system, i.e. 
when plugging £ = L/2 into (1^ . the entanglement en¬ 
tropy S{L/2) is expected to be a linear function of InL 
with a slope pc. The numerical data at the maximum of 
S{L/2) is in good agreement with this prediction as shown 
on hgure [H This confirms the divergence of the correla¬ 
tion length as the lattice size and therefore, the occurrence 
of a phase transition. Because of the magnetic and elec¬ 
tric fields coupled to the boundaries of the system during 


(26) 
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= 1/4 


e = 1/2 


£ = 1 



h 

£ = 2 




£ = 4 




-)<- 

-e- 


L = 8 
L = 10 
L = 12 
L = 16 
L = 20 


Fig. 6. Variance of disorder fluctuations of polarization. The different graphs correspond to different values of e and the different 
curves to different lattice sizes L. 


the numerical computations, the CFT predictions for the 
slope of S{L/2) with InL do not apply. 


4 Autocorrelation functions 

As discussed in the introduction, the average connected 
autocorrelation functions A{t) display three different be¬ 
haviors according to the values of the parameters S and 
e, i.e. the position in the phase diagram. On the critical 
lines, a slow relaxation ([3]) depending on the logarithm of t 
is expected. In the Griffiths phases, rare regions induce an 
algebraic decay m of the autocorrelation functions, with 
an exponent 1/z. Finally, away from the Griffiths phases, 
a more usual exponential decay is recovered. 

The algebraic decay of autocorrelation functions in the 
Griffiths phases has been observed numerically in the case 


of the random quantum Ising chain by exploiting the map¬ 
ping onto a gas of free fermions [14]. The lattice sizes that 
we were able to reach with DMRG being much smaller, 
such an algebraic decay of the spin-spin or polarization- 
polarization autocorrelation functions could not be ob¬ 
served for the random Ashkin-Teller model. A purely alge¬ 
braic behavior is indeed expected to hold only in the large¬ 
time limit t ^ 1 and in the thermodynamic limit L ^ 1. A 
transient regime may be observed for small times t while, 
for large times t, the finite-size of the system may induce 
an exponential decay of autocorrelation functions. Usu¬ 
ally, one looks for an intermediate regime in the numerical 
data where the asymptotic behavior holds. No such inter¬ 
mediate regime could be found for both the spin-spin or 
polarization-polarization autocorrelation functions. This 
is particularly clear when plotting an effective exponent 

d^in versus t. For values of the transverse field h ex- 
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Fig. 7. Entanglement entropy for an equal partition of the system, i.e. I = Lj^. The different graphs correspond to different 
values of t and the different curves to different lattice sizes L. 


pected to be in the Griffiths phases, two non-algebraic 
regimes, where the effective exponent varies with t, are 
observed at short and large times. But in between, no 
plateau corresponding to a purely algebraic decay could 
be distinguished. 


To fit our numerical data, we used an extended expres¬ 
sion of the one proposed by Rieger et al. for autocorrela¬ 
tion functions in a Griffiths phase m- The assumptions 
are the same: in the paramagnetic phase, the probability 
of an ordered region of linear size i scales as p{i) ^ 
and its tunneling time is t{£) ^ ^. In a finite system of 

width L, the linear size of rare regions is bounded by h so 


the average autocorrelation function reads 

A(t) = e-‘/^ = f 

Jo 

f — cjc' pt 

= -— / -^e-^dv 

= (27) 

where v = cr'/c = z is the dynamical exponent, and 

7 ( 0 , 3 :) is the incomplete gamma function. In the limit of 
large time t and lattice size L, one recovers the prediction 
^{0 = obtained in the saddle-point approxi¬ 

mation. 

The numerical estimates of the connected autocorrela¬ 
tion functions were fitted with the 4-parameter non-linear 
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Fig. 9. Spin-spin (left) and polarization-polarization (right) autocorrelation functions of the random quantum Ashkin-Teller 
chain versus time t. The different graphs correspond to different values of the transverse field h. The continuons fines correspond 
to a fit, either with the ansatz (1281) or with an exponential. In the legend, the number between parenthesis after the transverse 
field h is the estimate oiljz given by the fit. When the best fit is the exponential, exp is indicated instead. 



Fig. 8. Scaling of the maximum of the entanglement entropy 
S(Z//2) for an equal partition of the system, i.e. I = T/2, with 
the logarithm of the lattice size. The different symbols corre¬ 
spond to different values of e and the straight lines are linear 
fit of the data. 

ansatz 

A{t) = ait““=|7(a2,a3t) - 7(a2,a4<)l- (28) 

The bounds 0 < 02 < 1 were imposed during the fit¬ 
ting procedure. The quality of the fit was quantified us¬ 
ing the mean-square deviation ■ Because the bound¬ 
aries of the Griffiths phases are not known with a good 
accuracy, the data were also fitted with an exponential 
A{t) = Spin-spin and polarization-polarization 

autocorrelation functions are plotted respectively on fig¬ 
ures [9] for various transverse fields /i at e = 4. The contin¬ 
uous lines correspond to the best fit, Eq. (051) or exponen¬ 
tial, i.e. the one with the smallest mean square deviation 


X^. The inverse 1/z of the dynamical exponent is indicated 
in the legend when (1281) is the best fit, while exp indicates 
a fit with an exponential. As can be seen on the figures, 
the data are nicely reproduced by an exponential decay for 
large and small transverse fields h. Close to the transition 
point, the best fit is obtained with (l28ll . which means that 
the corresponding range of transverse fields is in a Griffiths 
phase. As expected, when e < 1, these phases are centered 
around h = 1 and their boundaries are similar for spin- 
spin and polarization-polarization autocorrelation func¬ 
tions. For e = 4, the Griffiths phase is shifted to smaller 
values of the transverse field for spin-spin autocorrelation 
functions and to larger ones for polarization-polarization 
autocorrelations. For e = 2, the shift is seen only for the 
polarization-polarization autocorrelation functions. It was 
also the case for the peak of the autocorrelation time (fig¬ 
ure El. At the boundaries of the Griffiths phases, the data 
is not well fitted, neither by an exponential form nor by 
the ansatz (EH- When the best fit is obtained with the 
ansatz (1281) . the dynamical exponent takes a value z = 1, 
i.e. saturating the imposed bound z < 1. The deviation 
between the fit and the numerical data is clearly visible 
on figures El The transverse fields for which such a devi¬ 
ation occurs, are probably in a region of cross-over where 
the autocorrelation functions display a more complex be¬ 
havior. 

On figures [TO] and [TTl the inverse of the dynamical 
exponents z is plotted versus the transverse field h. As 
conjectured in Ref. EH, the dynamical exponents display 
a peak centered at the corresponding critical point, i.e. 
at the magnetic transition for the dynamical exponent of 
spin-spin autocorrelation functions and electric transition 
for polarization-polarization autocorrelations. As already 
observed for other peaked quantities, the two transitions 
occur at the same control parameters for e < 1 and are 
separated for e > 1. Note that the maxima of the dy¬ 
namical exponents are found at the locations of those of 
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Fig. 10. Inverse of the dynamical exponent 2 estimated by a fit of the spin-spin antocorrelation functions with Eq. (I28II and 
plotted versus the transverse field h. The different curves correspond to different lattice sizes and the different graphs to different 
values of e. 


the autocorrelation times and of the variance of disorder 
fluctuations. Between these two transitions lines, there is 
therefore a double Griffiths phase, i.e. a disordered Grif¬ 
fiths phase in the magnetic sector and an ordered one in 
the electric sector where both dynamical exponents z are 
larger than 1. However, as seen on figures fTOl and fTTl these 
Griffiths phase are not infinite but have a finite extension, 
because of the binary distribution of the couplings Ji and 
Ki. For e = 4, it is observed that the magnetic and electric 
Griffiths phases still overlap. Nevertheless, it will proba¬ 
bly not be the case anymore for larger values of e. 

In the random Ising chain, the dynamical exponent was 
shown to behave as z ~ 1/2|d| in the Griffiths phase [JT]. 
A similar behavior seems to be also reasonable in the 
case of the random Ashkin-Teller chain, as can be seen 
on figures [in] and dn The boundaries = — In and 


(5_ = — In h- of the Griffiths phase were first estimated 
respectively as the first and last points with a dynami¬ 
cal exponent z > 1. The critical point is assumed to be 
located at 6c = (d+ -I- 6-)/2. The two dashed lines plot¬ 
ted on figures [ini and [TT] simply correspond to straight 
lines l/z{5) = {5 — 5c)/{6+ — 5c) for 5 G [dc;<5-i-] and 
l/z{5) = (5 — 5c)/{5- — 5c) for 5 G [d_;dc]. The slope 
is not equal to two, as in the Ising model, but is in the 
range 1 — 1.5. As the lattice size is increased, the numerical 
data seem to accumulate on these straight lines, at least 
at the boundaries of the Griffiths phase. In the neighbor¬ 
hood of the critical point, much larger lattice sizes would 
be necessary to test this linear behavior of 1/z. In the 
case of L = 20, the dynamical exponent seems to be over¬ 
estimated for h < he- Again, this may be explained by 
the under-sampling already observed with disorder fluc¬ 
tuations of magnetization and polarization. 
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e = 1/4 


e = 1/2 


£ = 1 



Fig. 11. Inverse of the dynamical exponent 2 estimated by a ht of the polarization-polarization autocorrelation function with 
Eq. (I28II and plotted versus h. The different curves correspond to different lattice sizes and the different graphs to different 
values of e. 


5 Conclusions 

The random quantum Ashkin-Teller chain has been stud¬ 
ied by means of time-dependent Density Matrix Renor¬ 
malization Group. The average over all possible disorder 
configurations was performed for L < 16. For L — 20, 
a partial average is observed to induce an under-sampling 
of disorder fluctuations of magnetization and polarization. 
Such partial averages are commonly used in the literature 
in the study of random systems. Our data show that they 
should be considered with great care, especially in the 
quantum case. 

The analysis of integrated autocorrelation times and of 
the variance of disorder fluctuations leads to a phase dia¬ 
gram qualitatively in agreement with the one conjectured 
by Hrahsheh et al. on the basis of SDRG m- However, 


finite-size effects are large, especially for entanglement en¬ 
tropy, and our lattice sizes are too small to allow for an 
accurate extrapolation in the thermodynamic limit. The 
coincidence of the maxima of disorder fluctuations with 
the critical lines confirms that the phase transition is gov¬ 
erned by disorder fluctuations, and not by quantum fluc¬ 
tuations. Nevertheless, the divergence of the entanglement 
entropy as the logarithm of the lattice size is recovered, as 
in pure quantum chains. In the regime e > I, the existence 
of a double Griffiths phase is confirmed. Using an origi¬ 
nal method to take into account finite-size effects, the two 
dynamical exponents, associated to the algebraic decay 
of spin-spin and polarization-polarization autocorrelation 
functions respectively, could be computed. They display 
the expected behavior in a Griffiths phase: a peak cen¬ 
tered at the magnetic or electric transition. Furthemore, 
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it seems reasonnable to assume that they diverge in the 
thermodynamic limit as z{6) ~ 

It is our pleasure to gratefully thank Cecile Monthus for dis¬ 
cussions and for having pointing out some useful references on 
the topic. 
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